% MLE estimation of parameters
% modified for RAND round 1 revision on June 25, 2022
% RUN this on RIVANNA on Interactive mode.
% Runs in within 1 hours.
clear
clc
deflator = 1; % 1 for basic deflation 2 for ppi of oil drilling

rng(5758934)

thisFile = mfilename('fullpath');
thisDir = fileparts(thisFile);
cd(fileparts(fileparts(thisDir)));

data0 = open('data/processed/OilData.mat');
dataQ = table2array(data0.OilData.Production)*30/1000;     % dataQ is in million barrels per month
dataQ_names = data0.OilData.Production.Properties.VariableNames;

addpath('code/2_estimation/f4_estimating_parameters_aux_files');

if deflator==1
    dataP = data0.OilData.Price; % Deflator (2019 Q4 =100): EXPLAIN IT IN THE PAPER
else
    dataP = data0.OilData.Price_oil_drilling_ppi_deflated;
end


T = size(dataQ,1); I = size(dataQ,2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% De-trending
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% detrending the data (x must be row vector)
tempfun = @(x,tdata) -exp(repmat(x(1:I),T,1)-x(2*I+1)*tdata)...
    + exp(repmat(x(I+1:2*I),T,1));
tdata = repmat((1:T)',1,I);

tempfun1 = @(x) sum(sum((dataQ - tempfun(x,tdata)).^2));
options = optimoptions('ga','MaxGenerations',250*(2*I+1));
tau_est = ga(tempfun1,2*I+1,[],[],[],[],zeros(1,2*I+1),[9^999*ones(1,2*I) 1],[],[],options);

%% IMPORTANT: hereafter we use detrended dataQ!
dataQ = dataQ + exp(repmat(tau_est(1:I),T,1)-tau_est(2*I+1)*tdata) ;

% compute means and std. deviations
dataQ_means = mean(dataQ);
dataQ_std = std(dataQ);
table0 = [dataQ_means' dataQ_std'];

%%% clustering in ng groups
% starting point of k-means algorithm (marked in graph)
%start_obs = [2 7 17 19]; ng = length(start_obs);

%[groupid, centroids] = kmeans(table0,ng,'MaxIter',10000,'Start',table0(start_obs,:));

% already did the Kmeans... then ussed KMEANS + geograpphy to get the
% following

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% ESTIMATION %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 6 groups
ng = 6;

% new classification that combines K-means and Geography
groupid = [4;... %Angola
    1;... %UAE
    2;... %Canada
    2;... %China
    4;... %Algeria
    3;... %Ecudaor
    4;... %Egypt
    2;... %UK
    1;... %Iran
    1;... %Iraq
    1;... %Kuwait
    4;... %Libya
    3;... %Mexico
    4;... %Nigeria
    2;... %Norway
    1;... %Qatar
    5;... %Russia
    5;... %Saudi
    6;...  %USA
    3];  % Venezuela

data1 = [dataP dataQ];
conv_ma = (repmat(groupid,1,ng)==repmat(1:ng,I,1));

% started with this initial parameter
if deflator==1
    param_ini = ...
        [0.2/30    0.1/30   47.9592  513.2512...
        0.3942    0.3851    0.2459  0.3851    0.2459  3.3491...
        0.9347    0.6093    0.8824  0.6093    0.8824 10.9169...
        0.0015    0.0039   18.8001];
else
    param_ini = ...
        [0.0414814988858829	1.00132639260428e-07	110.654641966071	1000	11	8.25013061969832	9.50287982893483	11	3.53890836259071	3.80546584014374	5.27955352413246	3.85701584659068	3.76255260254769	4.16995634790022	4.49836593633120	7.44550173715047	1.00018103292796e-07	1.00000629458240e-07	38.008884190512];
end


%call the estimation functions
tempres0 = estimate_fun(data1,conv_ma,I,ng,param_ini); % runs fminsearch 4 times
tempres = estimate_fun_midaco(data1,conv_ma,I,ng,tempres0(1:end-2)); %uses fminsearch outcome for midaco search 4 times
param_est = tempres(1:end-1);
llfvalue = tempres(end);

%%


theta_mle      = [param_est(1), ... % beta
    param_est(2), ... % lambda
    param_est(3), ... % muU
    param_est(4), ... % sig2U
    param_est(5:ng+4),... %alv
    param_est(ng+5:2*ng+4), ... % betv
    param_est(2*ng+5), ... %a1-tilde
    param_est(2*ng+6), ... %a2-tilde
    param_est(2*ng+7), ... %bar-w
    param_est(end), ... %ulow
    ];

if deflator==1
    save('output/estimates/ESTIMATES.mat','theta_mle','-v7.3')
else
    save('output/estimates/ESTIMATES_oil_gas_drilling_deflator.mat','theta_mle','-v7.3')
end

